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\ ABSTRACT 

o 
o 

Rotochemical heating originates in a departure from beta equihbrium due to spin-down com- 
' pression in a rotating neutron star. The main consequence is that the star eventually arrives 

at a quasi-equilibrium state, in which the thermal photon luminosity depends only on the cur- 
rent value of the spin-down power, which is directly measurable. Only in millisecond pulsars 
the spin-down power remains high long enough for this state to be reached with a substantial 
luminosity. We report an extensive study of the effect of this heating mechanism on the thermal 
^ \ evolution of millisecond pulsars, developing a general formalism in the slow-rotation approxima- 

• tion of general relativity that takes the spatial structure of the star fully into account, and using 

a sample of realistic equations of state to solve the non-superfluid case numerically. We show that 
I nearly all observed millisecond pulsars are very likely to be in the quasi-equilibrium state. Our 

' predicted quasi-equilibrium temperatures for PSR J0437-4715 are only 20% lower than inferred 

, from observations. Accounting for superfluidity should increase the predicted value. 

Subject headings: stars: neutron — dense matter — relativity — stars: rotation — pulsars: 

o 



tin 

in 



general — pulsars: individual (PSR J0437-4715 — PSR J0108-1431) 



1. INTRODUCTION 

, Neutron star cooling is an important tool for the study of dense matter. By comparing cooling models 

5-H \ with observations of thermal emission from these objects, we can gain insight into the equation of state (EOS) 

of dense matter, the signatures of exotic particles, superfluid energy gaps, and magnetic field properties [for 
a recent review, see Yakovlev & Pethick (2004)]. A neutron star loses the thermal energy with which it was 
born initially through neutrino emission, and later through photon emission, the change between these two 
stages occurring about 10^ yr after birth. 

Several heating mechanisms may become important during late stages of the evolution of a neutron 
star, most of them related to the delayed adjustment to the progressively changing equilibrium state as 
the rotation slows down. These include the dissipation of rotational energy due to interactions between 
superfluid and normal components of the star (Shibazaki & Lamb 1989) and release of strain energy stored 
by the solid crust due to spin-down deformation (Cheng et al. 1992). 

Another of these mechanisms is rotochemical heating (Reisenegger 1995), which has its origin in devi- 
ations from beta equilibrium due to spin-down compression. As a neutron star spins down, the centrifugal 
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force acting on each fluid element diminishes, changing the local value of the pressure. Since the composi- 
tion of neutron star matter in beta equilibrium is a one-parameter function, this compression results in a 
displacement of the equilibrium concentration of each particle species. Reactions which change the chemical 
composition are in charge of driving the system to the new equilibrium configuration. But if the rate at 

which reactions do this task is slower than the change of the equilibrium concentrations due to spin-down 
compression, the system is permanently out of chemical equilibrium. This implies an excess of energy, which 
is dissipated by enhanced neutrino emission and heat generation. 

After its introduction, this heating mechanism was studied by several authors: Cheng & Dai (1996) ap- 
plied it to quark stars, Reisenegger (1997) made an order-of-magnitude estimate of the effects of superfluidity, 
and lida & Sato (1997) studied the heating due to compositional transitions in the crust due to spin-down 
compression. All studies agree that rotochemical heating is particularly important for old neutron stars with 
fast rotation and low magnetic fields, features that are characteristic of millisecond pulsars (MSPs). 

The most striking prediction associated with rotochemical heating is that, if the spin-down timescale is 
substantially longer than any other timescale involved (with the likely exception of magnetic field decay), 
the star arrives at a quasi- equilibrium state, in which the temperature depends only on the current, slowly 
changing value of (the product of the angular velocity and its time derivative), proportional to the 
spin-down power, and not on the star's previous history (Reisenegger 1995). This provides a simple way 
to constrain the physics involved in theoretical models, once the spin parameters and observed surface 
temperature of a MSP are known. 

Although the qualitative behavior of the thermal evolution of neutron stars with rotochemical heating 
is known, all previous studies made only order-of-magnitude estimates, ignoring the spatial structure of 
the star and thus not oflfering reliable predictions to be compared with observations. Also unknown is the 
dependence of the thermal evolution on EOS and stellar mass. 

In this work, we calculate the thermal evolution of MSPs with rotochemical heating, taking the structure 
of the star fully into account in the frame of general relativity, and using realistic EOSs of dense matter. In 
order to do so, we develop a general formalism to treat the evolution of the temperature and departures of 
chemical equilibrium. As a first approach, we have chosen the simplest possible core composition, namely 
neutrons, protons, electrons, and muons (npe/i matter), to see if we can explain observations with it, before 
invoking exotic particles. The main difficulty of treating the spatial structure of the star is the lack of an 
expression for the spin-down compression in the frame of general relativity. We develop such an expression 
with the aid of the slow-rotation approximation of Hartle (1967). The astrophysical situation to be modelled 
is the stage after the accretion-driven spin-up of the pulsar has finished. 

Although it is generally believed that crusts and cores of neutron stars have superfluid components (e.g., 

Yakovlev & Pethick 2004), in this work we ignore superfluidity. Our plan is to make an extensive study of 
rotochemical heating in non-superfluid stars and compare results with observations, in order to assess whether 
invoking superfluidity (whose parameters are currently very uncertain) is necessary to explain observations. 

The structure of this paper is the following: In §2, we develop our general formalism, as well as the 
expression for the spin-down compression rate. In §3, we detail the input necessary for numerical calculations 
(EOS, reactions, heat capacities, etc.). In §4, we describe our results, comparing our predictions with the 
recent detection of likely thermal ultraviolet emission from the nearest MSP, PSR J0437-4715 (Kargaltsev, 
Pavlov, & Romani 2004), and the lack of optical emission from the even closer "classical" pulsar PSR J0108- 
1413 (Mignani, Manchester, & Pavlov 2003). We also list some other MSPs for which the rotochemical 
heating luminosity might become detectable. A short summary of our main conclusions is given in §5. 
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2. THEORETICAL FRAMEWORK 
2.1. Basic Equations 

It is conventional in neutron star cooling calculations to divide the star into two regions: a nearly 
isothermal interior, which ranges from the center out to density pb ~ 10^° g cm~^, and a thin envelope, 
which ranges from pb to the surface and where a strong temperature gradient exists (Gudmundsson, Pethick, 
& Epstein 1983). Since wc are modelling the thermal evolution of a MSP long after accretion has stopped, 
it is safe to assume that thermal relaxation from an initial non-uniform internal temperature profile has 
already occurred, so that the redshifted internal temperature, 

To, = r(r)e*M, (1) 

is uniform (see, e.g., Glendenning 1997). Here, gu = — e^* is the time component of the metric of a non- 
rotating reference star, of which r is the radial spherical coordinate. Although spherical symmetry is broken 
for a rotating star, we describe the latter as a perturbation to the corresponding non-rotating star (with the 
same total baryon number). We establish a Lagrangian correspondence between the surfaces of constant r 
in the non-rotating star with the constant-pressure surfaces of its rotating counterpart, on which all local 
thermodynamic quantities will be shown to be constant (see § 2.2). The evolution of the internal temperature 
is given by the thermal balance equation (Thorne 1977), which for an isothermal interior is given by 

f^ = ^[L^-L^-L^], (2) 

where C is the total heat capacity of the star, is the total power released by heating mechanisms, is 
the total neutrino luminosity, and is the photon luminosity. These quantities are calculated as 

c = E / '^^c^.- (3) 

IS = JdVQHe^^, (4) 

= J dVQ.e^'^, and (5) 

L~ = AnaR^Ty- =4naRl{T,^)\ (6) 

respectively, where dV = Airr'^ y/g^dr is the proper volume element, cv.i is the spcicific heat (heat capacity 
per unit volume) of each particle species, Qi/ is the total neutrino emissivity contributed by reactions, Qh 
is the total heating rate per unit volume, a is the Stefan-Boltzmann constant, R is the stellar coordinate 
radius, $s = ^(^)) ^oo = Re~^' is the effective radius as measured from infinity, and the redshifted 
effective temperature. The surface temperature T, is obtained from the internal temperature by assuming 
an envelope model (Gudmmidsson et al. 1983; Potekhin et al. 1997). 

Since the neutrino emissivity and heating rate are modified when the neutron star is out of chemical 
equilibrium (Haensel 1992; Reisenegger 1995), the evolution of the temperature depends on how strongly it 
departs from the beta equilibrium state. For npe/i matter, this departure can be quantified by the chemical 
imbalances (Haensel 1992): 



L 



Vnpe = Sfln - ^jlp - 6 fig, 



(7) 
(8) 
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where 6fii = /Zj — ^^'^ is the deviation from the cquihbrium chemical potential of species i, at a given pressure. 
Since diffusion timescales are short compared with the evolutionary timescales to be considered (Reisenegger 
1997), we assume uniform redshifted chemical potential deviations throughout the core: 

<5/x~=5Mi(r)e*M. (9) 

To obtain the time evolution of the chemical imbalances, we start by writing down the chemical potential 
of each particle species as a function of the number density of all particles: /ij = /ii({nj}). We assume 
small departures from chemical equilibrium, which can be expressed by requiring that <C In this 
approximation, the departures from the equilibrium particle number densities 6ni = rii — n^' are related to 
the by 

where the partial derivatives are evaluated at the beta equilibrium state. To eliminate the effect of particle 
diffusion between different regions of the star, we integrate equation (10) over regions where free particles 
exist. After integration, wc obtain the deviation from the equilibrium number of particles 6Ni as a function 
of the redshifted chemical potential lags: 

dNi = Y,BijStif, (11) 

j 

with 

B., = / dV^e-", (12) 

J core 

where we have used equation (9). Since the Bij do not depend on time, we invert and take the time derivative 
of equation (11), obtaining the evolution of the (5/xf°: 

6fir = T.iB-'),^5Nj. (13) 



The rate of change of SNi is given by 
where 



N, 



SNi = Ni- TVf , (14) 
Vi= [ dVe^Y^ATi (15) 

J core „, 



is the change in the total number of particles of species i due to reactions (Thornc 1977). Here, AF^ is the 
net creation rate of particles of species i per unit volume due to reaction a. It can be checked that the iVj 
satisfy baryon number and charge conservation. The quantity N^'^ is the change in the equilibrium value of 
the total number of particles of species i, which depends on the spin-down compression rate. We defer its 
calculation to the next subsection. The evolution of the redshifted chemical imbalances r]°° = 77(r)e* follows 
from equations (7), (8), and (9): 

VZe = - - ^f'T (16) 

Cp^ = ^Ajr - - (5A^ (17) 



Equations (2), (16), and (17) give a complete description of the thermal evolution of a neutron star with 
rotochemical heating and npeii composition, given an expression for N^'^ . 
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2.2. Lagrangian Spin-down Compression Rate 

For slow enough rotation frequencies, the deviation of the star from the non-rotating configuration can 
be treated as a small perturbation, which we describe in terms of a Lagrangian formalism. Using the fact that 
the total baryon number ^ of a star is conserved, we can describe its interior in terms of surfaces of constant 
pressure P enclosing a fixed number of baryons N (e.g., the surface P = encloses A baryons inside it). As 
the stellar rotation rate changes, these surfaces will readjust their shapes, with a corresponding change in 
P, but keeping A'' constant. Since there is a one-to-one relation between the enclosed number of baryons A'' 
and the radial coordinate of the non- rotating configuration r (for fixed A), we use them interchangeably to 
identify a given surface, relating them by dN = Airr'^ y/g^ndr . 

The underlying assumption is that, on a surface of constant pressure, the other thermodynamical quan- 
tities are also constant. This is straightforward for neutron-star matter in beta equilibrium (described by a 
barotropic equation of state), but not so obvious for departures from this state (which are not necessarily 
barotropic, as the pressure-density relation will also depend on the local particle abundances). We show in 
Appendix A that, in a uniformly rotating, perfect-fiuid star in hydrostatic equilibrium, surfaces of constant 
pressure coincide with those of constant energy density and gravitational redshift, regardless of whether 
the equation of state is barotropic. -'^ This, together with diffusive equilibrium (imiform redshifted chemical 
potentials everywhere in the stellar interior), ensures that the chemical composition on each surface is also 
constant. 

To quantify the spin-down compression rate and thus N^'' , we need an expression for the change in 
pressure with rotation frequency on each surface enclosing a constant N, keeping A constant. Hartle (1967) 
developed a perturbative approach in the frame of general relativity, relating quantities on rotating and 
non-rotating stars with the same central energy density pc (but, therefore, difi'erent A), by expanding the 
metric of an axially symmetric non-rotating star in even powers of the angular velocity fl. This perturbative 
approach is valid for stars rotating at frequencies much smaller than the Kepler frequency CIk, at which 
mass shedding from the equator occurs. For realistic equations of state, this limiting frequency can be well 
approximated by means of empirical formulae (Lasota, Haensel, & Abramowicz 1996). The corresponding 
Kepler periods, Pk, are lower than the shortest MSP periods measured to date [see, e.g.. The ATNF Pulsar 
Catalogue^ (Hobbs et al. 2004)]. Our task is to express our spin-down compression rate at constant A in 
terms of Hartle's results, which assume constant pc- 

If the baryon number enclosed by a surface of constant pressure P in a rotating star is A' = N{A, P, f2^), 
some partial-derivative manipulation yields 



dP \ (dN/dn^) 



P.A 



1 S(—\ -(^\ {dA/dn 



2\ 

Pc 



(18) 



where the subscripts denote quantities which are kept constant in carrying out the differentiation. We note 

that {dA/d^l'^)p^ and {dA/dpc)Q2 are just {dN/dfl'^)p^p^ and {dN/dpc)p^u^, respectively, evaluated at the 
surface of the star (P = 0). Also, since we are perturbing around the non-rotating configuration, we can 
evaluate all the derivatives at = 0. 



^The validity of the assumption of a perfect fluid is also discussed in Appendix A. 
^http://www.atnf.csiro.au/research/pulsar/psrcat/expert.html 
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Hartle (1967) fomid that the lowest-order change in a given stellar property with rotation is linear in 
the squared angular frequency , for a star with fixed central density. We can thus write 



dN 



N{p„P,n^)-N{p,,P,0) 
Q2 



(19) 



and consider it independent of Q^. To find the number of baryons enclosed by a surface of constant pressure, 
we recall the baryon number conservation law (Misner & Sharp 1964), and write 



N 



[ i-gy 

Jv 



(20) 



where g is the determinant of the metric, u* is the time component of the fluid's four-velocity, n is the 
rest frame baryon number density, and V is the spatial volume enclosed by the constant pressure surface 
(throughout this section, we use G = c = 1). For the non-rotating case, we have 



N{pc,P,0) 



47ry2e^(''-«)n(pc,y)dy, 



(21) 



where gr 



„2A 



(1 — 2M/r) ^ is the radial component of the non-rotating metric. For the rotating case. 



one has to express the perturbed metric of Hartle (1967) in the non-rotating coordinate system. Keeping 
terms to order O^, we get 



{-9) 



-1/2 



e*+^r2 sin( 



1 



2M 



.e-f($ + A) + ^+2fc+|i' 
dr r or 



1- 



i-*^^ L 1 ^2 -2/1-2 -2<I> 

t— h+ -r sm Ouj e 

dr 2 



(22) 
(23) 



where m, h, and k are metric-perturbing functions which are proportional to fi^ and can be separated into 
spherical and quadrupolar parts, u> is the difference between the star's rotation frequency and the frame- 
dragging frequency, and ^ is a Lagrangian displacement that relates surfaces of constant p in rotating and 
non-rotating stars of the same pc (Hartle 1967). The spherical polar angle is 9, and the mass enclosed inside 
the radius r is M. Replacing equations (22) and (23) in (20), keeping terms to order O^, taking the angular 
integral, and integrating by parts in r, we get^ 

N{p,,P,n^) = r^'"''^7ryV(^-^)n(p„2/) 1 



mo 



1 2-2 -2* 

-y u> e 
3" 



^o(Pc:?y) / dn\ 
n{Pc, y) \dyj ^ 



y-2M 

dy + 47rr^e'^(''='''^^o(Pc, r)n{pc, r). 



(24) 



The functions w, mo, are obtained by solving first and second-order differential equations in the coordinate 
r, as explained in Hartle (1967), and which need a value for the rotation rate of the star fl and a non-rotating 
relativistic stellar structure as input. The zero subscript on functions indicates their spherical part (the 
quadrupolar part vanishes when performing the angular integral, and the function k has no spherical part). 

The remaining derivatives in equation (18) can be obtained from the non-rotating configuration. We 

first have 



ON 

'dP 



02, A 



Anr'^e'^n 
dP/dr ' 



(25) 



^Hartle (1967) and Hartle & Thorne (1968) give an expression for the total number of baryons of the rotating star in which 
the last term of equation (24) is not present, which is correct when integrating up to the surface, but not at intermediate layers. 
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whcrc dP/ dr can be calculated analytically from the relativistic stellar structure equations (Oppenheimer & 
Volkoff 1939). The other derivative has to be obtained by differentiating equation (21): 

opcj \opc 



+ / 47ry' 



dPcJy y \dpc 



dy. (26) 



For numerical calculations, it is convenient to write the derivative in the first term of the right-hand side as 

dr\ _ idP/dp,)r ^27) 



in order to calculate the three derivatives with constant r in the same way. 

The resulting spin-down compression rate, equation (18), is plotted in Figure 1, normalized by pressure 
and Kepler frequency. The left panel displays results obtained for different EOSs, with fixed central pressure, 
while the right panel shows results for different stellar models and fixed EOS. As expected, the expression 
is negative, since spin-down compression means increase in P with decrease of . Compression is strongest 
at the center of the star. It can be shown that the central value is given by 



dP \ Pclcipc + Pc) (dA/dn 



2\ 



J Pa 



J na'^^^ PcIT'c {dA/dpc)n^ 



(28) 



where 7 = d log P/ d log p is the adiabatic index, and the subscripts c mean central values. Since {dA/dpc)Q2 = 
for the maximum mass non-rotating configuration, the divergence of the result for M M^ax, as shown 
in the right panel of Figure 1, is easy to understand. The graphs confirm the order-of- magnitude expression 

which implies that the perturbation is small as long as f2 <C fix- 

With our expression for the spin-down compression rate at hand, we can write the change in the 
equilibrium number of particles of species i with time as 

■ r dV"^ f dP\ 

Nt" = 2nn / dN-^ —2 = 2nnin,i, m 

J core \CSi / N,A 

where the derivative dY,^ /dP is calculated in beta equilibrium. It is straightforward to check that the N^'^ 
satisfy baryon number and charge conservation. 



2.3. Effect of a Phase Transition 

The variables N and P are clearly continuous functions of r in the non-rotating star, regardless of the 
presence of a phase transition. In Appendix B, we show explicitly that the same is true for {dP/d^^)N.A- 
On the other hand, the equilibrium particle concentrations may have a discontinuous jump AY^'^ at the 
transition, which would contribute a Dirac delta function to dY^'^/dP. Therefore, there would be a finite 
contribution to N^'' from the transition, 



m / dP ^ 

dP\di^yN,A 



(31) 



- 8 - 



To order of magnitude, the fractional contribution of a phase transition to the total value of N^'^ (and 
to the total heating of the star, if the reactions are of the same kind everywhere) is the ratio of the jump 
AY^*' to the total variation in Y^'^ across the stellar core. In the present work, as described in the following 
sections, phase transitions are not important, as the only transition in one of our EOSs is fairly weak^. In 
other cases, such as a star with a core of deconfined quarks of fairly uniform relative abundances, the phase 
transition may dominate the heating. 

In order to maintain diffusive equilibrium, a finite number of particles per unit time, Nl'^\trans, must 
be brought to (or, if negative, away from) the phase transition much faster than the latter progresses due 
to the spin-down of the star. Therefore, a phase transition is the place where our assumption of diffusive 
equilibrium is most strongly put into question. However, in the case of interest to us, the low temperatures 
(implying long mean-free paths and therefore fast diffusion) and the extremely slow evolution should still 
validate this approximation. 

Glendenning (1992) has shown that, due to the presence of two separate conserved charges (baryon 
number and electric charge), phase transitions in neutron stars may not occur as a single discontinuity, 
but instead through a mixed phase, in which both phases coexist over a finite range in pressure (or stellar 
radius). In this case, any discontinuities would be considerably weakened, if not fully eliminated, and we do 
not expect any special effects (such as discussed above) associated with the phase transitions. 

3. INPUT FOR THE NON-SUPERFLUID CASE 

To implement the formalism of §2 numerically for the non-supcrfluid case, we need as input an EOS, 
neutrino emissivities, a specific heat, and an envelope model. We then have to calculate the stellar structure 
and the corresponding integrals over the star. 

3.1. Equations of State 

Since our formalism requires a detailed knowledge of the structure of the star, wc need an EOS to deter- 
mine the spatial dependence of the thermodynamical quantities by solving the relativistic stellar structure 
equations (Oppenheimer & Volkoff 1939). In addition to the usual variables given in most published EOSs 
(pressure, energy density, chemical composition, and adiabatic index, all evaluated in chemical equilibrium), 
we need to know the partial derivatives of equation (10) (which are easy to compute if we know the energy 
density or energy per baryon of interacting particles as a function of baryon number density and proton 
fraction) and the baryon effective masses for computing reaction rates and heat capacities. We also want 
to cover a broad range of EOSs, in order to assess the dependence of the quasi-equilibrium temperature on 
their properties. Given our requirements, we have used three sets of EOSs: 

1. the two most realistic EOSs from Akmal, Pandharipande, & Ravenhall (1998) (APR) for the core, 
supplemented with that of Pethick, Ravenhall & Lorentz (1995) and Haensel & Pichon (1994), for the 
inner and outer crust, respectively; 



*In our numerical integration, the discontinuity is automatically taken care of by evaluating the derivative dY^'^ /dP in each 
integration step as the ratio of the (not necessarily small) increment in Y'"' to the small increment in P. 
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2. five representative EOSs from Prakash, Ainsworth & Lattimer (1988) (BPAL) for the core, with the 
same EOSs for the inner and outer crust as for the previous set; and 

3. a non-interacting Fermi gas (see, e.g., Shapiro & Teukolsky 1983) for the whole star. 

In Table 1 we list, for each EOS, the mass M, central energy density pc, coordinate radius R, and effective 
radius as seen from infinity i?oo of the maximum-mass non-rotating configuration. 

Since in the APR set the thermodynamical quantities are obtained from an analytical fit to tabulated 
many-body calculations (Akmal et al. 1998), care has to be taken when interpreting results obtained for 
densities higher than the highest tabulated value, since extrapolation can lead to significant errors. This is 
the reason why we do not list the maximum-mass configuration for the A18 + 6v EOS. Another feature of 
this set is that the A18 + Sv + UIX* EOS becomes non-causal at densities greater than 2 x lO^^gcm"^, 
which is the central density of a star of 2.14M0, slightly below the maximum mass formally allowed for this 
EOS, M^ax = 2.19M0. However, this EOS is considered to be the state of the art among those derived from 
non-relativistic potentials, since it incorporates three- nucleon interactions and relativistic boost corrections 
(Akmal et al. 1998). For this reason, we will in what follows take it as our reference EOS. 

The EOSs of the BPAL set are labelled according to Prakash et al. (1997), who characterize them by the 

different magnitude of their nuclear incomprcssibility and density dependence of their symmetry energy. The 
values for the Fermi gas EOS listed in Table 1 correspond to a central density slightly below the appearance 
of S~ hyperons. 

The A18 + Sv + UIX* EOS has an additional feature: a phase transition associated with the appearance 
of a neutral pion condensate at a density 4 x 10^^ g cm~^. In order to connect the normal phase with the 
pion-condcnscd one, we assumed a Maxwell transition, which implies a discontinuity in all thermodynamical 
quantities but the temperature, pressure and neutron chemical potential, resulting in an energy-density jump 
of 6.6%. Akmal et al. (1998) also considered the possibility of a mixed-phase transition (e.g., Glendenning 
1992), concluding that (for a standard, IAMq neutron star) the shell containing the mixed phase is only ^ 40 
m thick. They also estimated the Debye screening length of electrons to be smaller than the characteristic 
size of structures (by a factor 6-12), which would favor the Maxwell transition as a better approximation to 
the real situation (Heiselberg, Pethick, & Staubo 1993). 

Also listed in Table 1 are the Kepler periods Pk — 27r/i}K corresponding to most of the EOSs, which are 
necessary in order to determine the extent to which the slow-rotation approximation holds for treating the 
spin-down compression. The APR and PAL values were computed with the empirical formiila of Lasota et al. 
(1996). This gives, for each EOS, the maximum allowed rotation frequency in terms of the mass and radius 
of the maximum-mass non-rotating configuration. The accuracy for non-causal EOSs is somewhat lowered 
from its optimal value of 1.5% (Lasota et al. 1996), which applies for A18 +5v + UIX*. The reason why 
the entry corresponding to A18 + Sv is empty is that its maximum mass configurations occur at a central 
density which requires substantial extrapolation, making the empirical formula an unreliable estimator of 
Pk- For the Fermi gas EOS, we adopt the exact value calculated for a pure neutron gas by Haensel, Salgado, 
& Bonazzola (1995), since the empirical formula of Lasota et al. (1996) does not work for this case. So far, 
the fastest known MSP is PSR B1937-)-21, with a period of 1.56 ms (Backer et al. 1982), which is about 
three times longer than Kepler periods for all realistic EOSs. 
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3.2. Neutrino Emissivity 

Since we are not taking superfluidity into account, direct and modified Urea reactions are the dominant 
neutrino emission mechanisms. We have therefore neglected additional neutrino emission processes in the 
core and the crust of the star, and ignored possible modifications in Urea rates due to the neutral pion 
condensate of the APR set (see, e.g., Khodel et al. 2004). 

Away from beta equilibrium, neutrino emissivities and net reaction rates per unit volume of Urea-type 
reactions in non-superfluid matter are given by (Haensel 1992) 

Q„(n,T,ry„) = g^«(n,T)F. (^) (32) 

AT„{n,T,dr,a) = ^Ql'{n,T)H. , (33) 

respectively, where is the neutrino emissivity in equilibrium due to reaction a, rja is the chemical 
imbalance affected by reaction a, T is the local temperature, n is the baryon number density, k is Boltzmann's 

constant, and and ff* are dimensionlcss control functions which depend on whether the reaction is of 
direct of modified Urea type. These functions arc given by (Rciscnegger 1995) 

^ , , , 1071x2 315a;'' 2lx^ 

= '^^^^^^^^ ^''^ 

, , 714a; 420a;3 42a;5 

"^^^"^ = 457^^+ 457^+ 457^ ^"^'^ 

^ , , . 22020a;2 5670a;4 420a;*5 ^x^ ,„„, 

^^^"^^ = ^+ 11513;? + Ti5i3;^^ + n5i3;?^+m ^^^^ 

14680X 7560a;3 840x^ 24a;^ 
^^""^ " 11513;?^ + 11513^ + 11513;? + 11^^ ^^^^ 

where the subscripts D and M mean direct and modified Urea, respectively. A finite r\a of either sign 
enhances neutrino emission due to the even nature of the functions i^* . The amount of energy released by 
each reaction of type a is r]a (Reisenegger 1995), thus the total energy dissipation rate per unit volume is 



Qff = ^Ar„77„, (38) 

with the signs defined by 

Ar = Va^b-Vb^a (39) 

■q = S^{A) - 5ii{B), (40) 

where Ta^b is the rate per unit volume of the reaction that transforms the particle set A to the set B. This 
sign convention ensures that reaction rates are enhanced in the direction which restores equilibrium, since 
the functions are odd. 

The equilibrium emissivities of both direct and modified Urea reactions can be written as 

Q^«(n,T) = 5„(n)r«, (41) 

where Sa is a slowly varying function of n, g = 6 for direct Urea, and q = S for modified Urea (e.g., Yakovlev 
et al. 2001). This decoupling of the temperature from the spatial part introduces a great simplification. To 
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obtain the non-equilibrium neutrino luminosities and net reaction rates, we first define 

(42) 

r2e^5a(n)e*(2-9)dr, (43) 

where we have used equations (1), (5) and (9), and where i'^eg the equilibrium neutrino luminosity due 
to reaction a. Combining equations (32), (33), (42) and (43), we obtain 
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I 



= L^F4^^)Tl, (44) 
AT^e^dV = ^L„H4U)T^-\ (45) 



We denote modified Urea reactions with electrons and muons by a = Me and a = Mjj,, respectively, in 
each case adding the contributions of the neutron and proton branches (Yakovlev et al. 2001). Direct Urea 
with electrons or muons is denoted by a = De and D/z, respectively. From the APR set, the only EOS 
which allows direct Urea is A18 + Sv + UIX*, exceeding the threshold for direct Urea with electrons at 
PDe = 1-59 X 10^^ g cm~"^, which is the central density of a star of 2 Mq. Direct Urea with muons lies in the 
non-causal regime. With the exception of BPAL31, all the EOSs of the PAL set allow electron direct Urea. 
Muon direct Urea is forbidden for both BPAL21 and BPAL31. 

Combining equations (4), (9), (33), (38), (42), and (45), we obtain the heating luminosity corresponding 
to the reaction a 

L^,„=L^iMia)Tl. (46) 

We can combine the heating and cooling contributions of each reaction into a single expression, which will 
prove to be useful for physical insight. If we define 

M,(ea) = ^aH4^^) - F,(^„), (47) 

we can write the difference of equations (44) and (46) as 

L^.c,-L^ = LM^c)Tl. (48) 

Thus, the relative size of rj^ and /cTqo will determine, through the functions M*, whether there is net cooling 
or heating due to the neutrino emitting reactions. Since the functions M* are even, they do not depend on 
the sign of the We plot the functions M*(^) for the direct and modified Urea case in Figure 2. With 

^ = 0, the constant term in K,, gives the conventional cooling case. As |^| grows, the cooling is enhanced by 
the additional neutrino emission due to the chemical imbalance, reaching a maximum cooling at |^| ~ 3.5. 
For larger values of |^|, the heating becomes important, completely balancing the cooling for |^| ~ 5.5. For 
very large values of |^|, the heating dominates, growing as ^ and ~ in the direct and modified Urea 
case, respectively. In this limit, a fixed fraction of the energy released is emitted as neutrinos, 3/8 and 1/2 
for modified and direct Urea, respectively. The remainder stays behind, heating the star. 



3.3. Envelope 

Since we want to simulate the evolution of a MSP which has finished accretion, we use the fully accreted 
envelope model of Potekhin et al. (1997) to calculate the effective surface temperature 

= lO^V (l.81e-*''Too,8)' '' K^ (49) 
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where gi4 is the surface gravity in units of 10"'^'' cm s~^, Too,8 is the internal temperature in units of 10^ 
K, and = $(rb). The photon luminosity then follows from equation (6). An accreted envelope is more 
heat transparent than a pure iron one for > 10^ K, leading to faster cooling, whereas for Tg < 10^ K the 
opposite happens (Potekhin et al. 1997). 

3.4. Heat Capacity 

To calculate the integral in equation (3), we use the specific heat capacities at constant volume cv,i for 
degenerate, non-superfluid fermions (see, e.g., Levenfish & Yakovlev 1994). In analogy with the neutrino 
emissivity, we can separate the temperature dependence from the spatial one. From equation (3), we can 
define 

C=-^ = ^Y.j ^^r^e''m*{n)pFi{n)e-^dr, (50) 

i 

where ppi and m* are the Fermi momentum and effective mass of species i, respectively. The latter can be 
obtained analytically for the APR and PAL EOSs (see, e.g., Page et al. 2004), while for the noninteracting 
Fermi gas we use m* = /ii/c^, where /i^ is the chemical potential of species i. The sum is performed over all 
particle species and the integral is done over the region where each particle species exists. We take only free 
particles into account, neglecting the heat capacity due to the lattice of ions in the crust. 

3.5. Crustal Processes 

In this work, we arc neglecting processes in the crust which might also modify the total number of 
neutrons, protons, or electrons. We may ask here if this choice is a good approximation for computing 
the evolution of the chemical imbalances, since the short diffusion timescale limit means that the chemical 
equilibrium c;very where in the star can be restored by reactions occurring wherever these are fastest. lida & 
Sato (1997) calculated the heating of a neutron star due to spin-down compression in the crust. They found 
that neutron absorptions in the inner crust are the dominant non-equilibrium process, producing a very small 
heating rate H{t) 5 x 10~^£^, where E is the spin-down power. Although they did not consider neutron 
diffusion between different layers in the star, and modelled the spin-down compression as a one-dimensional 
process, we take their result to suggest a very small contribution to the total heating. Since the inner crust 
accounts for about 2% of the total baryons in the star and the neutrino emissivities of the core overwhelm 
crustal emission processes, we integrated equations (12) only over the core and neglected processes in the 
crust, under the assumption that the error introduced would be small. Our results confirm our assumption, 
since for the more realistic EOSs we obtained heating rates due to core processes much bigger than lida & 
Sato's results (see equation 69). We have to warn that neglecting the crust may not be a good approximation 
if superfluidity is included, since non-equilibrium processes in the crust can even overtake suppressed Urea 
reactions in the core. It has also been suggested by Gusakov et al. (2004) that direct Urea reactions could 
occur at the crust-core interface. We have not explored this possibility. However, we are aware that even a 
small amount of direct Urea reactions in the inner crust would have a significant effect on the evolution of 
chemical imbalances. 
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3.6. Thermal Evolution Equations 

Finally, we write the coupled equations for the time evolution of the temperature and the chemical 
imbalances explicitly for the non-superfluid case with npen composition, taking direct or modified Urea 
reactions into account. For the internal temperature, from equations (2), (3), (4), (5), (6), (48), and (50), 
we get 



T = 



1 

+ 



{MD{^npe)LDe + M D{^npy.)L D ^) 

For the chemical disequilibria, using equations (11), (14), (15), (16), (17), and (30), we obtain 



(51) 



fi°° 



npe 



InpL 



where we have defined 



(^L DeHDi^npe)T^ + LMeHMi^npe)T^^ 
LDfiHD{^npfj,)T^ + LMnHM{inpii)Tl^ + 2Wnpe^^^^, 
{LDeHD{^npe)T^ + LMeHM{^npe)T^^ 

- (^LD,iHDi^npn)T^ + LMnHMi^npi^)T^^ + 2W„p^nf2, 



Bn 



Bn 



Bnn-^pp BjipBpn 



Be 



(52) 
(53) 

(54) 

(55) 
(56) 



and 



^^npe — i^^npe ^np^^Q.e ~t~ Zjipl^^ p^ 
^^npfj. ~ {^npfi ^np^^Cl,{j. H" Z^plo^^p. 



(57) 
(58) 



Typical values for the different constants are plotted in Figure 3. Since VlVl < 0, equations (52) and (53) 
have a positive term, proportional to the spin-down power, which arises from the change in the equilibrium 
concentrations of each particle species with spin-down and makes the chemical imbalances grow. The re- 
maining terms, opposite in sign to ry^g and ?7^^, account for the effect of reactions trying to restore beta 
equilibrium. Equation (51) has the photon luminosity as a negative term, and the net contribution of neu- 
trino reactions, equation (48), which may be positive or negative depending on the absolute value of ^npe and 
Cnpiu- ^'^^ numerical calculations, the evolution of flQ was computed assuming magnetic dipole braking with 
no field decay, relating magnetic field, rotation period, and period derivative by the conventional formula 
B 2:^3.2 X 10^^ V PP G, where P is measured in seconds. 
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4. RESULTS AND DISCUSSION 

4.1. Thermal Evolution 

A typical solution of equations (51), (52), and (53) is shown in the left panel of Figure 4. First, the 
temperature starts falling as in the conventional cooling case, while the chemical imbalances grow due to 
spin-down compression. At some point, the ratios ^npe and ^„p^ (right panel of Figure 4) are big enough, so 
that the net contributions of neutrino reactions, equation (48), are positive, and their sum exceeds L'^ . The 
temperature then starts rising. The chemical imbalances continue to grow, more slowly than the temperature 
and hence reducing ^npe and ^„pp- Finally, the star can arrive at a quasi- equilibrium, state, where the rate 
at which spin-down modifies equilibrium concentrations is the same as the rate at which reactions drive the 
system towards the equilibrium configuration, with heating and cooling balancing each other (Reisenegger 
1995). The subsequent thermal evolution can then be approximated by the simultaneous solution of the 
equations Too = and fi^^ = j)^^ = 0. We discuss this quasi- equilibrium solution in detail in §4.2. 

The pure modified Urea case {Ljje = Ld/j, = 0) is the simplest to analyze. Figure 5 shows the thermal 
evolution of the same star as in Figure 4, this time with different initial conditions of internal temperature 
(left) and chemical imbalances (right). The spin-down parameters were chosen so that the arrival at the 
quasi-equilibrium state is clearly visible, the short-dashed lines being the quasi-equilibrium solutions. We 
note that, in both cases, the value of the temperature at the quasi-equilibrium state and the time required 
to arrive do not depend on initial conditions. 

If direct Urea reactions are open, solutions change noticeably. In Figure 6 we show the thermal evolution 
of three different stellar models built with the same EOS, with identical initial conditions and spin-down 
parameters. The IAMq star is the same as that of Figure 5, with only modified Urea reactions in operation. 
The 2Mq star has electron direct Urea operating, with the corresponding muon reaction forbidden. We see 
that the minimum temperature is higher and occurs at a much earlier time, after which the system reaches a 
"metastable" quasi-equilibrium state, corresponding to partial equilibration between Too and r]^^. However, 
^rap/i continues to grow, leading the system to a final quasi-equilibrium state, determined by muon modified 
Urea. The 2.17M0 star has direct Urea with electrons and muons. This time the final quasi-equilibrium 
state occurs even earlier than the "metastable" quasi-equilibrium of the 2Mq star, with the temperature at 
about the same level. 



4.2. The Quasi-Equilibrium State 

When only modified Urea reactions operate, it is possible to solve analytically for the quasi-equilibrium 
values of the photon luminosity L^^^^ and chemical imbalances ??^e,eg and rj^^^^g, as function of stellar 
model and current value of iltl. What makes these approximations possible is the high value of ^npe and 
£,nptj. at quasi-equilibrium. From the right panel of Figure 4, we see that, at quasi-equilibrium {t 10^'^ 
yr), ^npe ^ ^npij, ^ 200. This enables us to ignore all but the greatest power in the functions Hm and Mm 
(equations 36, 37, and 47): 

24 

^ Tl513;^5^'^^^^' 
15 

Mm{x) ~ 11513^8 ^' = (60) 



The error in this approximation, for x 200, is less than one part in 100. 
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To solve for i^eg' Too = V?^e — = and use equations (59) and (60) to rewrite (51), (52), 

and (53) as 

^Me(Cpe,e«)' + ^M4<^,e,)' = k^L^/CM (61) 
ZnpeLMeiV^pe^eq) + ^npLn ^liv'^p^^.q) = — W„pefifi (62) 

ZnpLMe{V?tpe,eq) + ^npiJ.LMfi{ri?tp^,eq) = -^Wnpn^l (63) 

Solving for rj^pg^eq V?^piJ,,eq ill equations (62) and (63), we get 

1 /7 

VZe,eq = k { (64) 




Cp^,e, = M7rT^I W (65) 

where we have used equations (57) and (58) and the fact that /n,p = In,e + Iii,n- Both chemical imbalances 
at quasi-equilibrium are positive, since f2f2 and are negative. We finally replace equations (64) and (65) 
in (61), to obtain the photon luminosity: 



^l,eq -'^M 



/8 \'/' //« \'/'' 
Lmb I \ Lmii J 



\Q.tl\^l\ (66) 



An interesting feature of equation (66) is that it does not depend on envelope model. The reason is that 
the limit ry°° fcToo implies that reactions are completely determined by chemical imbalances, eliminating 
any dependence on internal temperature. The energy released is divided in fixed fractions between neutrino 

emission (1 — \Cm/Ch] = 3/8) and heating (Cm/Ch = 5/8), the latter being radiated entirely as photons. 
Note also that the ratio Ti^e,eq/v^ij.,eq depends only on EOS and stellar model, and is of order unity, thus 
we will make reference to either of them as i]'^ . 

For the entire range of EOSs and stellar models, we can write 

, ^ 1030-31 (^) ' erg S-, (67) 

y ^ms J 

where P-20 is the period derivative measured in units of IQ-^", and Pms the period in ms. Values for different 
stellar models built with the AI8 + 5v + UIX* EOS are shown in Figure 3. Using equation (6), we can 
rewrite equation (67) in terms of the effective temperature: 




r~ :..(2-3)xlOM-3^^ K. (68) 



This time, the uncertainty duo to EOS and stellar model is smaller than 25%. We can also write the ratio 
of quasi-equilibrium photon luminosity (and thus total heating rate) to the spin-down power E 

L°° / P \ 

^'^^ - ■ (0.3-3) X 10-M . (69) 



E ' [PS. 



ms 
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It is also possible to calculate ^£^eq for the quasi-equilibrium state. Its spin-down power dependence, 

^e,e<i oc (70) 

where a = 2.42 is the exponent of the envelope model of Potekhin et al. (1997) (see §3.3), shows that 
increasing |f2f2| reduces ^e, decreasing the accuracy of the analytical approximation, although very slowly. 



4.3. Effect of Hyperons 

Although we did not include hyperons in our calculations, wc may ask how their presence coiild modify 
the thermal evolution. The first particles to appear after muons are probably the S~ and A'' hyperons, 
which require the introduction of two additional chemical imbalances: 

llnnSp = 2^„ — Ms — Mp) (71) 
VnA = Hn-HK- (72) 

Once hyperons are present, the following non-leptonic reactions are open (e.g., Langer & Cameron 1969): 

n + n ^ + p, (73) 
n-|-A° ^ (74) 

Reaction (73) proceeds via weak interactions, since it does not conserve strangeness, while (74) proceeds via 
strong interactions. Both reactions have timcscalcs several orders of magnitude shorter than beta processes 
(Langer & Cameron 1969). Therefore, imbalances (71) and (72) remain small compared to rinpe and r]nptj,, 
and reactions (73) and (74) contribute negligibly to the total heat generation. Moreover, since chemical 
imbalances associated to direct or modified Urea reactions with hyperons are linear combinations of TjnnHpi 
VnA, Vnpe, and rinpn, the latter two will determine the heat generation through both nucleon and hyperon 
reactions. 

In order to assess the importance of including the Urea processes involving hyperons in addition to the 
nucleon processes, consider, as an example, the S~ direct Urea reactions 

T,^^n + e + i>e, n + e^'E~ + D(., (75) 

whose associated chemical imbalance is rj-^ne = Vnpe — Vnn'Sp ~ Vnpe- Their net effect on equations (52) 
and (53) is to enhance the electron direct Urea rate according to L'^^ = Ljj^ + L^ne^ reducing the chemical 
imbalance and thus the stellar surface temperature. This correction is small, since direct Urea reactions 
with hyperons are at least a factor 5 weaker than their nucleon analogs (Prakash et al. 1992). For the pure 
modified Urea case, it can be checked from equation (66) that the correction to the surface temperature 
introduced by adding several reactions involving hyperons is a factor [i„/(I/„ + Lh)]^^"^^, where L„ and 
are the nucleon and hyperon Urea luminosities, respectively. Thus, the corrections due to hyperons in a 
purely modified Urea or purely direct Urea scenario are fairly negligible. 

The only important effect of hyperons in the context of rotochemical heating is that, in their presence, 

conditions for direct Urea processes may become more easily satisfied, as is the case with reactions involving 
A° (Prakash et al. 1992). The steep increase of the proton concentration with the appearance of S~ hyperons 
(e.g., Glendenning 1997) can also cause the condition for nucleon direct Urea to be satisfied at lower densities 
than if hyperons are excluded from the models. 
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4.4. Conditions for Arrival at the Quasi-Equilibrium State 

We may now ask which of the known MSPs are Hkely to be in the quasi-equilibrium state. To answer 
this question, we need to know how long it takes to arrive at the quasi-cquilibrium state, and if real pulsars 
are older than this time. Since the negative "equilibration" terms due to reactions in equations (52) and 
(53) are important only for imbalances near their equilibrium values (recall the limit ^ 1), we can assume 
that the initial evolution of jj^ is due only to Qfl: 

•nf « 2Wi^tl. (76) 

Integrating over time, we get 

r^T{t)^m[nl-n\t)\. (77) 

For the true solution to have reached the quasi-equilibrium state, this approximate solution must exceed 
?7£^g(Of2). Using equations (64) or (65), together with equation (77), we may express this condition as an 
upper limit on the initial spin period 

P 

^/f+A 

with 

1 (2k^In,tV" {ntlfl^ 



Po < = Pt^ (78) 



m V chU j ■ ^^^^ 

(Here and in the rest of this subsection, we assume that only modified Urea processes are active.) Given 
the current spin parameters of a pulsar, condition (78) gives the highest initial period it can have had, so 
that enough time has elapsed for it to have reached the quasi-equilibrium state. For MSPs, the constant A 
is generally small. For a 1.4Mq star built with the A18 + Sv + UIX* EOS, we get 

A^om{p.2oP^y/'. (80) 

Values for other stellar masses are shown in Figure 3. Prom equation (76), we can estimate a characteristic 
timescale for equilibration as 

„oo / p3 \ 6/7 

= 1.6 X 10^ ( 4^ ) yr. (81) 



' 2Winn \P-2o 

Reisenegger (1995) calculated the thermal evolution of neutron stars with rotochemical heating for 
different values of the magnetic field, using the magnetic dipole braking model with no field decay. His 
conclusion was that a necessary condition for the quasi-equilibrium solution to be a good approximation 

to the exact one is that the spin-down timescale has to be longer than that for the growth of chemical 
imbalances. To quantify this, we estimate the spin-down timescale as the characteristic age, 

P ^ 

and take the ratio of Teq to it, getting 

^=A. (83) 

Tsd 

Figure 7 shows the internal temperature of our conventional star after the rise from the minimum tempera- 
ture, for different initial values of A, with fixed initial T^q. Dotted lines are the quasi-equilibrium solutions 
at the instantaneous value of 00. The figure suggests that, for high values of A, exact and quasi-equilibrium 
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solutions depart from each other, whereas for small A, solutions remain "together". In fact, all solutions 
depart from each other, the departure being slower with time for smaller A. For more general spin-down 
laws with arbitrary braking index n = f20/f2^, it can be shown analytically that this is true for 1 < n < 13. 
The increase in log (T^^/T'^^), where T®^ is the exact solution and T** the quasi-equilibrium one, is roughly 
linear with time. For A = 0.25 at t = 0, and assuming a standard braking index n = 3, T'^'^ is ^ 5% 
lower than T'^^ at t = T^q, the difference growing 2% every 10^ yr. Rewriting equation (80) in terms of the 
magnetic field, 

A « 0.01(P^,B8)'/', (84) 

where Bg, is the surface magnetic field in units of 10^ G, we can easily understand the results of Reisenegger 
(1995), since increasing B with fixed Pq increases the initial value of A. Relating A at the present time and 
at t = for an arbitrary, but constant, braking index, 

. p X (13-n)/7 

A{P,P)={^—j A{Po,Po), (85) 

shows that the former is greater than the latter for all reasonable values of n. Thus, a small current value 
of A is a good indicator that the star remains close to the quasi-equilibrium state, as long as it satisfies 

Po < pr- 

Since constraints on the initial periods Pq"^ of a few MSPs have been obtained from the cooling ages 
of their white dwarf companions (Hansen & Phinney 1998), we can check which objects are likely at the 
quasi-equilibrium state. In Table 2, further discussed in §4.6, we show the range in P™'^ for some MSPs 
obtained by Hansen & Phinney (1998). We made no distinction between the n = 2 and n = 3 case, and 
adopted the widest possible range in Pg"^. In Table 2 we also show the value of A and the upper limit on 
initial spin period Pq^ required for the MSP to be currently in the quasi-equilibrium state. Among MSPs 
with constraints on initial periods, PSR J1G12+5307 is the only one not satisfying the conditions. The small 
values of A imply that any initial periods substantially shorter than the current ones will leave the MSPs in 
equilibrium by the present time. 



4.5. Predictions for PSR J0437-4715 and PSR J0108-1431 

Recently, Kargaltscv ct al. (2004) measured ultraviolet emission from PSR J0437-4715, the most prob- 
able explanation being that it corresponds to thermal radiation. Since this object satisfies P,^'' < Pg"^ with 
small A, we assume that it is already at the quasi-equilibrium state. In Figure 8, we plot the blackbody 
fit of Kargaltsev et al. (2004) as dashed lines, which correspond to 68% and 90% confidence contours. We 
overplot our values of T^^^ for the spin parameters of PSR J0437-4715 (van Straten et al. 2001), as function 
of i?oo, for the EOSs listed in Table 1. The bold lines show the range in Roo corresponding to the mass 
constraint of van Straten et al. (2001), Mpsr = 1.58±O.18M0. The complete mass range for each EOS goes 
from IMq to 0.95Mjnax- The upper limit was chosen to avoid the divergence in the spin-down compression 
rate for masses near M^ax (see Figure 1), which makes T^^^ also diverge due to /q^, (see equation 66). 

The best agreement with observations is reached if only modified Urea reactions are present. When 
direct Urea reactions open, there are two abrupt drops in T^^^ with increasing stellar mass, the first one 
(^ 10%) due to electron direct Urea, and the second (~ 50%) due to muon direct Urea. This can be seen in 
Figure 8 for the curves calculated with BPALll and A18 + Sv + UIX*, as examples of only electron direct 
Urea, and in the curves corresponding to BPAL32 and BPAL33, which have electron and muon direct Urea. 
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Perfect agreement could be obtained by allowing masses within ^ 1% of the maximum-mass non-rotating 
configuration, for the EOSs with only modified Urea. However, we consider this scenario very unlikely. 

Using the EOSs which have only modified Urea reactions within the allowed mass for PSR J0437-4715, 
we can constrain the predicted effective temperature to the narrow range T^^^ = (6.9 — 7.9) x 10^ K, about 
20% lower than the blackbody fit of Kargaltsev et al. (2004). There are three possible reasons why we are 
not matching their results: 

1. We are not taking superfluidity into account. This would reduce Urea reaction rates, lengthening the 
equilibration timescale and raising the quasi-equilibrium temperature (Reisenegger 1997). 

2. We are neglecting other heating mechanisms (some of them directly related to superfluidity), which 
could further raise the temperature at any stage in the thermal evolution. Nonetheless, in MSPs, all 
proposed mechanisms are less important than rotochemical heating (Schaab et al. 1999). 

3. The thermal spectrum could deviate from a blackbody, in the same way as the isolated neutron star 
RX J1856-3754, which has a well-determined blackbody X-ray spectrum that underpredicts the optical 
flux (Walter fc Matthews 1997), indicating a more complex spectral shape of its thermal emission. 

Kargaltsev and collaborators stress the fact that PSR J0437-4715, despite its much larger spin-down 
age, has a higher surface temperature than the upper limit for the younger, "classical" pulsar J0108-1431, 
Tg°° < 8.8 X 10^ K, inferred from the optical non-detection by Mignani et al. (2003). In this regard, we 
note that the spin-down power (c>c f2f2) of the latter pulsar is 680 times lower than that of J0437, making 
rotochemical heating substantially less important. Its equilibration timescale, according to equation (81), is 
2 X 10^^ yr, longer than the age of the Universe and certainly much longer than the spin-down age of the 
pulsar. Thus, it is not expected to be even close to reaching its quasi-equilibrium state (although it is old 
enough to have lost its initial heat content). Its actual temperature should therefore be substantially smaller 
than the predicted quasi-equilibrium surface temperature, which is already (680)^^^ ~ 6 times lower than 
that of J0437, well below the observational upper limit. 



4.6. Predictions for other millisecond pulsars 

In Table 2, we provide predictions of T^°^^ (assuming modified Urea reactions and no superfluidity) for 
several MSPs. Since the high-energy part of their expected thermal spectrum is highly absorbed by neutral 
hydrogen in the interstellar medium, we chose to order them by decreasing predicted quasi-equilibrium flux 
in the low-energy (Rayleigh-Jeans) regime, 

FRJ,e<i (X (86) 

where d is the distance to the object. Each Fjij^eq entry in Table 2 is scaled by the value for PSR J0437-4715. 
The first group of MSPs after PSR J0437-4715 are nearby, single MSPs, for which Faj^eq is expected to be 
larger than for the binary MSPs, for which the initial period has been estimated from their white dwarf 
companion (see §4.4). We consider the former as the primary targets for observations analogous to those of 
Kargaltsev et al. (2004). However, even for PSR J0437-4715, the nearest and brightest MSP known to date, 
the detection of thermal emission is rather difficult, so the observation of other objects of Table 2 will be a 
real challenge. The second group of MSPs are those with estimates for Pq"^. However, their low fluxes make 
the detection of their thermal emission impossible for current instruments. 
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Optical detections are even more daunting than those in the ultraviolet. The extrapolated Rayleigh- 
Jeans spectrum of PSR J0437-4715 gives magnitudes U = 26.8, B = 28.1, and V = 28.5, which are completely 
overwhelmed by the white dwarf companion. The most promising single MSPs have about three times lower 
expected fluxes at Earth, i. e., are another 1.2 magnitudes fainter, beyond the reach of current telescopes. 

5. CONCLUSIONS 

In this work, we have made an extensive study of the effect of rotochemical heating in non-superfluid 
millisecond pulsars. We have set up a general formalism for the thermal evolution in the framework of 
general relativity, which takes the stellar structure fully into account and can be extended to include both 
superfluidity and exotic particles. A key ingredient in this formalism is the spin-down compression rate 
based on the slow-rotation approximation of Hartle (1967). 

The main consequence of rotochemical heating is the arrival at a quasi-cquilibrium state, in which the 
effective temperature of a neutron star depends only on the current value of the spin-down power. We 
argue that most of the known MSPs are very likely in the quasi-equilibrium state. If only modified Urea 
reactions are allowed, the quasi-equilibrium bolometric photon luminosity in the non-superfluid case can be 
well approximated by 

, ^ 10^0-31 (^] erg s-\ (87) 

independent of the neutron star envelope model. 

The influence of EOS and stellar model on the quasi-equilibrium state is very weak, the only significant 
factor being the occurrence of direct Urea reactions. If they are open for both muons and electrons, quasi- 
equilibrium temperatures are low, as in the conventional cooling case. If they are open only for electrons, 
the system arrives at a "metastable" quasi-equilibrium, after which it proceeds as if only modified Urea 
reactions with muons were present. The highest temperatures are reached when all direct Urea reactions are 
forbidden. 

Even our highest predicted quasi-equilibrium effective temperatures are lower than the blackbody fit of 
Kargaltsev et al. (2004) to the UV emission of PSR J0437-4715 by about 20%. The inclusion of superfluidity 
will likely raise our predicted temperatures, being the subject of future work. 

We thank G. Pavlov and O. Kargaltsev for letting us know about their work in advance of publication 
and for kindly providing the data for Figure 8. The authors are also grateful to M. Taghizadeh, M. van 
Kerkwijk, R. Mignani, D. Page, M. Catelan, M. A. Diaz, C. Dib, and P. Jofre for discussions that beneflted 

the present paper, and an anonymous referee for thoughtful comments that improved its final version. This 
work made extensive use of NASA's Astrophysics Data System Service, and received financial support from 
FONDECYT through grant # 1020840. 

A. Coincidence of Surfaces of Constant Pressure and Energy Density in General Relativity 

In this Appendix, we show that, in uniformly rotating, rclativistic stars in hydrostatic equilibrium, 
composed of a perfect fluid, the surfaces of constant pressure and those of constant energy density coincide. 
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as in the Newtonian case, even if the equation of state is not barotropic. In addition, the gravitational 
redshift is constant on the same surfaces. This result, together with the diffusive equilibrium condition, 
allows us to describe the spin-down compression in terms of Lagrangian changes of thermodynamic variables 
on isobaric surfaces enclosing a fixed baryon number, since all thermodynamical quantities are constant 
on each isobar. As a by-product, we also provide an alternative derivation of the equation of hydrostatic 
equilibrium for uniformly rotating, relativistic stars. 

The assumption a perfect (i.e., non-dissipative) fluid, whose stress-energy tensor depends only on the 
energy density p and (isotropic) pressure P, is of course not strictly true in the astrophysical situation of 
interest in this paper, in which entropy is generated by non-equilibrium weak interactions, particle diffusion, 
and heat conduction, and lost (from the star) through the emission of neutrinos and photons. However, the 
total energy dissipated (and eventually lost) along the star's lifetime is much smaller than that associated 
with the mass, random motions, and interactions of its particles, as well as its rotation, and the time scale 
of the release of this energy is many orders of magnitude longer than the dynamical time and the rotation 
period of the star. Therefore, the contribution of these dissipative processes to the energy and momentum 
fluxes is extremely small and can be neglected in the stress-energy tensor. This is in line with the usual 
approach (also used in neutron star cooling calculations and in the rest of this paper) of neglecting thermal 
effects when calculating the structure of the star, which is regarded as a fixed background on which thermal 
processes take place. 

The metric of a stationary, axially symmetric system can be written in the form (see, e.g., Hartle 1967): 

ds^ = gttdt^ + Qrrdr'^ + geedO'^ + g4,4,d(f)'^ + 2gt4,dtd(f), (Al) 

where the metric coefficients are functions of r and 9 only. The components of the 4-velocity of any fluid 
element of the uniformly rotating star are = = 0, and u'^ = Om*, with 

«* = (- [gtt + 2ngt4, + n^g4,4] ) . (A2) 

We note that 1/u* is the gravitational redshift factor, which corrects the energy of a freely moving particle, 
as it is measured by observers co-moving with fluid elements on different surfaces inside the star.^ 

The energy-momentum conservation equation can then be written as 

T«^.^ = {p + P)u''.y + {p + P)u''u^.^ + P,^g''^ = 0, (A3) 

where T"^ are the components of the stress-energy tensor of a perfect fluid. Applying to this equation a 
projection operator orthogonal to the four- velocity, paf} = gai3 + Wa^/j (Schutz 1985), we get 

(5a/3 + u^up){p + P)u^V'l^u^ + P^ = 0. (A4) 



Since the only non- vanishing terms are those with 7 and 5 taking the values {t, (/)}, the relevant connection 
coefficients are 

0, if a = t, (j); 



'^'si — { _ln°'i^r' o ^^^^ 



Hence, equation (A3) yields 

-\<yP + P)u'u^95-y,0 = 0. (A6) 



^Thc locally measured energy is E = —u^^Pfj. = u^(—pt — ^P^,), where are the covariant components of the particle's 
4^momentum. Due to the symmetries of the metric, the expression within parenthesis is conserved along the particle's world 
line, and is its redshifted energy. 
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Using equation (A2), we arrive at the equation of hydrostatic equilibrium 

P/3-(p + P)(lnu*),^ = 0, (A7) 

which is in agreement with the result of Cook, Shapiro & Teukolsky (1992), for the special case of rigid 
rotation. This shows that the redshift factor I/m* is also constant on the surfaces of constant P. Thus, we 
can write 

= -{p + P), (A8) 



which implies that p is also constant on the same surfaces. 



B. Continuity of the Compression Rate across Phase Transitions 

In this Appendix, we show that each term of the spin-down compression rate, equation (18), is continuous 
across a phase transition, even if the energy density and baryon number density change discontinuously. We 
use units with G = c = 1. 

The discontinuous quantities in the derivative {dN/dP)Q2^j^, equation (25), are n and dP/dr. The Gibbs 
free energy per unit volume for npefi matter in beta equilibrium at T = is 

P + p = M„n. (Bl) 

Since, at a phase transition, the pressure P and neutron chemical potential /i„ are continuous, the fractional 
size of the number density jump can be written as 

^ = (B2) 

n P+p ^ ' 

where An is the positive jump in number density. From the rclativistic stellar structure equations, we have 

dP , ^.M + Attt^P 



The fractional size of the corresponding jump across the phase transition is 

AdP/dr _ Ap 
dP/dr ^ p + P' 



(B4) 



Thus, from equations (B2) and (B4), we see that the quantity n/(dP/dr) is continuous across a phase 
transition, and so is {dN/dP)Q2 y^. 

Regarding the derivative {dN/di}^)p^sot equation (19), since the jump in n is integrated in equation (21), 
N{pc, P, 0) is continuous. For N{pc, P, Cl^), the jump in n due to the last term of equation (24) is cancelled 
by the contribution from dn/dr in the integral. Using 

— = —An5{r — rt) + continuous part, (B5) 

dr 

where 5 is Dirac's delta function and rt = r{pc,Pt), with Pt the pressure at which the phase transition 
occurs, we can write 

dn 

- / 4wy'^e^^^^^o{y)—dy = 4nr^e^^''*^^o{rt)AnQ{r-rt)+contmuouspaxt, (B6) 

Jo 
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where 9 is the step function. Equation (B6) is equal to minus the jump (increasing r) of the last term in 

equation (24). 

Finally, for the derivative (dN / dpc) p_q2 , the contribution from {dn/dpc)r to the integral in equation (26) 
cancels out the discontinuity due to n in the first term. Using 

AnS{pc — Pc) + continuous part 

S{i — rt) + continuous part, (B7) 

r=rt 

where P{p*,r) = Pt, and ri{pc) the radial coordinate at which {dn/dpc)r = 0, we can write 

Jo \9pcjy \ri-rt\ 

which is the negative of the jump (increasing r) of the first term in equation (26), since the sign of {dr/ deo) p 
is also (rj - n)/|rj - n]. 

The right panel of Figure 1 shows the spin-down compression rate inside the core of several stellar models 

built with the A18 + 5v + UIX* EOS, which has a phase transition. The occurrence of the transition can 
be seen as a (fairly harmless) discontinuity in the derivative of each curve. 




An 



dr 
dpc 




Qir - n) + cont. part, (B8) 
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Fig. 1. — Left: Spin-down compression rate for four different EOSs, normalized by pressure and Kepler fre- 
quency. Stellar models have fixed central pressure Pq = 4.5 x 10'^^ dyn cm~^. Right: Fractional compression 
for different stellar models built with the A18 + dv + UIX* EOS. Stars are labelled by their mass in units 
of the maximum-mass non-rotating configuration (2.I9M0). 
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Fig. 2. — Functions Md and Mm, corresponding to equation (47) for direct and modified Urea reactions, 
respectively. 
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Fig. 3. — The four higher panels show results for constants defined in section § 3.6, obtained from different 
stellar models built with the A18 + Sv + UIX* EOS. The two lower panels show the dependence on stellar 
model of the quasi-cquilibrium photon luminosity (equation 66) and the constant A, assuming Pms = -P-20 = 
1 and using the A18 + Sv + UIX* EOS. For the latter, the two values are obtained by taking £ = npe and 
£ = npn in equation (79), respectively. 
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Fig. 4. — Evolution of the internal temperature and chemical imbalances (left) and the ratios £, = 77°°/ {kT^) 
(right), for a IAMq star calculated with the A18 + 5v + UIX* EOS, with initial conditions Too = 10^ K 
and null chemical imbalances at f = 0, and spin-down parameters B = 10^ G and Pq = 1 ms. 



log t [yrl 



log t [yr] 



Fig. 5. — Evolution of the internal temperature for different initial conditions on temperature (left) and 
chemical imbalances (right). For both plots we set r]^^ = r]^^ = r}°° at t = 0. Fixed initial values are 
rf° = (left) and Too = 10* K (right). The short-dashed line is the quasi-equilibrium solution, obtained 
by solving = and yj^g = 77^^ = 0. The stellar model and spin-down parameters are the same as in 
Figure 4. 



Table 1. Maximum- mass non-rotating configuration and Kepler period for the equations of state used in 

this paper 



EOS 


^max 


Pc 


R 


Roc 


Pk" 




(Me) 


(IQis g cm-3) 


(km) 


(km) 


(ms) 


A18 + 6v 


1.55^ 


1.86 


9.81 


13.42 




A18 + 5v + UIX* 


2.19'' 


2.78 


9.97 


16.79 


0.51 


BPAL 11 


1.42 


4.45 


8.42 


11.86 


0.54 


BPAL 21 


1.70 


3.46 


9.33 


13.69 


0.56 


BPAL 31 


1.91 


2.86 


10.10 


15.18 


0.59 


BPAL 32 


1.95 


2.66 


10.56 


15.64 


0.63 


BPAL 33 


1.97 


2.53 


10.92 


15.94 


0.66 


Fermi gas 


0.62^^ 


1.10 


12.77 


13.80 


0.98 



^Corresponds to maximum value tabulated in Akmal ct al. (1998) 

''Stellar model lies in the non-causal regime of this EOS 

'^Corresponds to maximum mass before appearance of hypcrons 

''Calculated with empirical formula (sec text), except the last value, 
which was adopted from Haensel et al. (1995) 



6 

log t [yr] 



Fig. 6. — Evolution of surface temperature for different stellar models, calculated using the A18 + Sv + 
UIX* EOS, with fixed initial conditions. The mass of each configuration is labelled on each curve. The 2Mq 
star is slightly above the threshold for direct Urea with electrons, but below the threshold for muon direct 
Urea. The 2.17Mq star is near the maximum-mass non-rotating configuration, and has direct Urea with 
electrons and muons. The spin-down parameters are the same as in Figure 4. 
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Fig. 7. Exact and quasi-cquilibrium solutions for the internal temperature of the same star of Figure 4 
(solid and dashed lines, respectively), with fixed initial temperature, and different magnetic fields B, choosing 
the initial period Po in each case so as to satisfy Teq = 1.5 x lO'' yr, and therefore implying different initial 
values of A. 
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Fig. 8. — Quasi-equilibrium effective temperatures T^^^ obtained with different EOSs and stellar models, 
for the spin parameters of PSR J0437-4715. Dashed lines are the 68% and 90% confidence contours of the 
blackbody fit of Kargaltsev et al. (2004) to probable thermal emission from this pulsar. Bold lines indicate, 
for each EOS, the range corresponding to the mass constraint of van Straten et al. (2001) for PSR J0437- 
4715, MpsR = 1.58 ± 0.18Mq. Abrupt reductions in temperature with (increasing or decreasing) radius 
correspond to opening of direct Urea reactions. 
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Table 2. Predictions for MSPs likely to be observable or with estimates for initial spin period 



Object 


P 


pa 




s,e.q 




A 





r>wd 



Refs. 




(ms) 


(10-20) 


(kpc) 


(10^ K) 


{Fo437,eg) 




(ms) 


(ms) 




J0437-4715 


5 


76 


1.86 


0.14 


0.72 


1 


0.17 


5.33 


2.4-5.3 


1,2 


J1024-0719 


5 


16 


<1.85 


<0.25 


<0.79 


~ 0.35'^ 


<0.14 


>4.82 




3,4 


J2124-3358 


4 


93 


1.23 


0.27 


0.74 


0.28 


0.13 


4.64 




3,4 


J0030+0451 


4 


87 


<1.00 


0.32 


<0.70 


<0.19 


<0.12 


>4.60 




4,5 


J1744-1134 


4 


07 


0.71 


0.36 


0.74 


0.15 


0.09 


3.90 




6 


B1257+12 


6 


22 


4.26 


0.45 


0.86 


0.12 


0.22 


5.63 




3,4 


J0034-0534 


1 


88 


<0.67 


0.53 


<1.41 


<0.14 


<0.03 


>1.85 


<1.4 


2, 4,7 


J1012+5307 


5 


26 


1.34 


0.41 


0.71 


0.11 


0.14 


4.93 


>5.1 


2, 4,8 


B1855+09 


5 


36 


1.74 


0.90 


0.76 


0.03 


0.15 


5.00 


<4.6 


2, 3, 9 


J1713+0747 


4 


57 


0.81 


1.10 


0.70 


0.02 


0.11 


4.35 


2.2-3.1 


2, 10 


J1640+2224 


3 


16 


0.16 


1.15 


0.60 


0.01 


0.05 


3.08 


>1.6 


2,4, 11 


J2019+2425 


3 


93 


<0.70 


1.49 


<0.76 


<0.01 


<0.08 


>3.87 


0.9-3.9 


2, 4, 12 



'^Intrinsic spin period derivative, based on the latest available distance and proper motion. In cases 
where these quantities are not known well enough for a reliable proper-motion correction, the measured 
period derivative is given as an upper limit. 

•^Parallax distance when available, otherwise dispersion-measure distance based on the NE2001 Galac- 
tic electron density model (Cordes & Lazio 2002). In case of PSR J1024-0719, upper limit based on 
measured proper motion and condition of positive intrinsic period derivative. 

'^Only rough reference value, since calculated as ratio of two upper limits. 
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